{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Elastic 3D beam structures\n",
    "\n",
    "This tour explores the formulation of 3D beam-like elastic structures in FEniCS. One particularity of this tour is that the mesh topology is 1D (beams) but is embedded in a 3D ambient space. We will also show how to define a local frame containing the beam axis direction and two perpendicular directions for the definition of the cross-section geometrical properties. This tour is illustrated on a shell-like structure made of beams and is available in `.xdmf` format."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "The complete source files including mesh files can be obtained from  The complete source files including mesh files can be obtained from :download:`beams_3D.zip`.\n",
    "\n",
    "We also refer to the following tours for similar problems on beams or plates: :ref:`BeamBuckling` or :ref:`ReissnerMindlinQuads`\n",
    "\n",
    "The final deformed structure (under self-weight) should look like this:\n",
    "\n",
    ".. figure:: shell_deformed.png\n",
    "   :scale: 50%\n",
    "   :align: center"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variational formulation: beam kinematics and consititutive equations\n",
    "\n",
    "The variational formulation for 3D beams requires to distinguish between motion along the beam direction and along both perpendicular directions. These two directions often correspond to principal directions of inertia and distinguish between stong and weak bending directions. The user must therefore specify this local orientation frame throughout the whole structure. In this example, we will first compute the vector $t$ tangent to the beam direction, and two perpendicular directions $a_1$ and $a_2$. $a_1$ will always be perpendicular to $t$ and the vertical direction, while $a_2$ will be such that $(t, a_1, a_2)$ is a direct orthonormal frame.\n",
    "\n",
    " In the following, we will consider a Timoshenko beam model with St-Venant uniform torsion theory for the torsional part. The beam kinematics will then be described by a 3D displacement field $\\boldsymbol{u}$ and an independent 3D rotation field $\\boldsymbol{\\theta}$, that is 6 degrees of freedom at each point. $(u_t,u_1,u_2)$ will respectively denote the displacement components in the section local frame $(t, a_1, a_2)$. Similarly, $\\theta_1$ and $\\theta_2$ will denote the section rotations about the two section axis whereas $\\theta_t$ is the twist angle.\n",
    "\n",
    "The beam strains are given by:\n",
    "\n",
    "* the **normal strain**: $\\delta = \\dfrac{d u_t}{ds}$\n",
    "* the **bending curvatures**: $\\chi_1 = \\dfrac{d\\theta_1}{ds}$ and $\\chi_2 = \\dfrac{d\\theta_2}{ds}$\n",
    "* the **shear strains**: $\\gamma_1 = \\dfrac{du_1}{ds}-\\theta_2$ and\n",
    "    $\\gamma_2 = \\dfrac{du_2}{ds}+\\theta_1$. \n",
    "* the **torsional strain**: $\\omega = \\dfrac{d\\theta_t}{ds}$. \n",
    "\n",
    "where $\\dfrac{dv}{ds} = \\nabla v \\cdot t$ is the tangential gradient. Associated with these generalized strains are the corresponding generalized stresses:\n",
    "\n",
    "* the **normal force** $N$\n",
    "* the **bending moments** $M_1$ and $M_2$\n",
    "* the **shear forces** $Q_1$ and $Q_2$\n",
    "* the **torsional moment** $M_T$\n",
    "\n",
    "The beam constitutive equations (assuming no normal force/bending moment coupling and that $a_1$ and $a_2$ are the principle axis of inertia) read as:\n",
    "$$\\begin{Bmatrix}\n",
    "N \\\\ Q_1 \\\\ Q_2 \\\\ M_T \\\\ M_1 \\\\ M_2 \n",
    "\\end{Bmatrix} = \\begin{bmatrix}\n",
    "ES & 0 & 0 & 0 & 0 & 0\\\\\n",
    "0 & GS_1 & 0 & 0 & 0 & 0\\\\\n",
    "0 & 0 & GS_2 & 0 & 0 & 0\\\\\n",
    "0 & 0 & 0 & GJ & 0 & 0\\\\\n",
    "0 & 0 & 0 & 0 & EI_1 & 0\\\\\n",
    "0 & 0 & 0 & 0 & 0 & EI_2\n",
    "\\end{bmatrix}\\begin{Bmatrix}\n",
    "\\delta \\\\ \\gamma_1 \\\\ \\gamma_2 \\\\ \\omega \\\\ \\chi_1 \\\\ \\chi_2 \n",
    "\\end{Bmatrix}$$\n",
    "\n",
    "where $S$ is the cross-section area, $E$ the material Young modulus and $G$ the shear modulus, $I_1=\\int_S x_2^2 dS$ (resp. $I_2=\\int_S x_1^2 dS$) are the bending second moment of inertia about axis $a_1$ (resp. $a_2$), $J$ is the torsion inertia, $S_1$ (resp. $S_2$) is the shear area in direction $a_1$ (resp. $a_2$).\n",
    "\n",
    "The 3D beam variational formulation finally reads as: Find $(\\boldsymbol{u},\\boldsymbol{\\theta})\\in V$ such that:\n",
    "$$\\int_S (N\\widehat{\\delta}+Q_1\\widehat{\\gamma}_1+Q_2\\widehat{\\gamma}_2+M_T\\widehat{\\omega}+M_1\\widehat{\\chi}_1+M_2\\widehat{\\chi}_2)dS = \\int \\boldsymbol{f}\\cdot\\widehat{\\boldsymbol{u}}dS \\quad \\forall (\\widehat{\\boldsymbol{u}},\\widehat{\\boldsymbol{\\theta}})\\in V$$\n",
    "\n",
    "where we considered only a distributed loading $\\boldsymbol{f}$ and where $\\widehat{\\delta},\\ldots,\\widehat{\\chi}_2$ are the generalized strains associated with test functions $(\\widehat{\\boldsymbol{u}},\\widehat{\\boldsymbol{\\theta}})$.\n",
    "\n",
    "## FEniCS implementation\n",
    "\n",
    "After first loading the shell mesh file in `.xdmf` format, we will use the UFL `Jacobian` function which computes the transformation Jacobian between a reference element (interval here) and the current element. For the present case, the Jacobian is of shape (3,1). Transforming it into a vector of unit length will give us the local tangent vector $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "from ufl import Jacobian, diag\n",
    "from plotting import plot\n",
    "\n",
    "mesh = Mesh()\n",
    "filename = \"shell.xdmf\"\n",
    "f = XDMFFile(filename)\n",
    "f.read(mesh)\n",
    "\n",
    "def tangent(mesh):\n",
    "    t = Jacobian(mesh)\n",
    "    return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(inner(t,t))\n",
    "\n",
    "t = tangent(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now compute the section local axis. As mentioned earlier, $a_1$ will be perpendicular to $t$ and the vertical direction $e_z=(0,0,1)$. After normalization, $a_2$ is built by taking the cross product between $t$ and $a_1$, $a_2$ will therefore belong to the plane made by $t$ and the vertical direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "ez = as_vector([0, 0, 1])\n",
    "a1 = cross(t, ez)\n",
    "a1 /= sqrt(dot(a1, a1))\n",
    "a2 = cross(t, a1)\n",
    "a2 /= sqrt(dot(a2, a2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define the material and geometrical constants which will be used in the constitutive relation. We consider the case of a rectangular cross-section of width $b$ and height $h$ in directions $a_1$ and $a_2$. The bending inertia will therefore be $I_1 = bh^3/12$ and $I_2=hb^3/12$. The torsional inertia is $J=\\beta hb^3$ with $\\beta\\approx 0.26$ for $h=3b$. Finally, the shear areas are approximated by $S_1=S_2=\\kappa S$ with $\\kappa=5/6$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "thick = Constant(0.3)\n",
    "width = thick/3\n",
    "E = Constant(70e3)\n",
    "nu = Constant(0.3)\n",
    "G = E/2/(1+nu)\n",
    "rho = Constant(2.7e-3)\n",
    "g = Constant(9.81)\n",
    "\n",
    "S = thick*width\n",
    "ES = E*S\n",
    "EI1 = E*width*thick**3/12\n",
    "EI2 = E*width**3*thick/12\n",
    "GJ = G*0.26*thick*width**3\n",
    "kappa = Constant(5./6.)\n",
    "GS1 = kappa*G*S\n",
    "GS2 = kappa*G*S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now consider a mixed $\\mathbb{P}_1/\\mathbb{P}_1$-Lagrange interpolation for the displacement and rotation fields. The variational form is built using a function `generalized_strains` giving the vector of six generalized strains as well as a function `generalized_stresses` which computes the dot product of the strains with the above-mentioned constitutive matrix (diagonal here). Note that since the 1D beams are embedded in an ambient 3D space, the gradient operator has shape (3,), we therefore define a tangential gradient operator `tgrad` by taking the dot product with the local tangent vector $t$.\n",
    "\n",
    "Finally, similarly to Reissner-Mindlin plates, shear-locking issues might arise in the thin beam limit. To avoid this, reduced integration is performed on the shear part $Q_1\\widehat{\\gamma}_1+Q_2\\widehat{\\gamma}_2$ of the variational form using a one-point rule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ue = VectorElement(\"CG\", mesh.ufl_cell(), 1, dim=3)\n",
    "W = FunctionSpace(mesh, Ue*Ue)\n",
    "\n",
    "u_ = TestFunction(W)\n",
    "du = TrialFunction(W)\n",
    "(w_, theta_) = split(u_)\n",
    "(dw, dtheta) = split(du)\n",
    "\n",
    "def tgrad(u):\n",
    "    return dot(grad(u), t)\n",
    "def generalized_strains(u):\n",
    "    (w, theta) = split(u)\n",
    "    return as_vector([tgrad(dot(w, t)),\n",
    "                      tgrad(dot(w, a1))-dot(theta, a2),\n",
    "                      tgrad(dot(w, a2))+dot(theta, a1),\n",
    "                      tgrad(dot(theta, t)),\n",
    "                      tgrad(dot(theta, a1)),\n",
    "                      tgrad(dot(theta, a2))])\n",
    "def generalized_stresses(u):\n",
    "    return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])), generalized_strains(u))\n",
    "\n",
    "Sig = generalized_stresses(du)\n",
    "Eps =  generalized_strains(u_)\n",
    "\n",
    "dx_shear = dx(scheme=\"default\",metadata={\"quadrature_scheme\":\"default\", \"quadrature_degree\": 1})\n",
    "k_form = sum([Sig[i]*Eps[i]*dx for i in [0, 3, 4, 5]]) + (Sig[1]*Eps[1]+Sig[2]*Eps[2])*dx_shear\n",
    "l_form = Constant(-rho*S*g)*w_[2]*dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clamped boundary conditions are considered at the bottom $z=0$ level and the linear problem is finally solved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "Matplotlib plotting backend only supports 2D mesh for scalar functions.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-24-6ced7cf47285>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mu\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mFunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mk_form\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0ml_form\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/Fenics/comet-fenics/examples/beams_3D/plotting.py\u001b[0m in \u001b[0;36mplot\u001b[0;34m(object, *args, **kwargs)\u001b[0m\n\u001b[1;32m    446\u001b[0m             cpp.log.info(\"Object cannot be plotted directly, projecting to \"\n\u001b[1;32m    447\u001b[0m                          \"piecewise linears.\")\n\u001b[0;32m--> 448\u001b[0;31m             \u001b[0mobject\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mproject\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobject\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmesh\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmesh\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    449\u001b[0m             \u001b[0mmesh\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mobject\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfunction_space\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmesh\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    450\u001b[0m             \u001b[0mobject\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mobject\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cpp_object\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Fenics/comet-fenics/examples/beams_3D/plotting.py\u001b[0m in \u001b[0;36m_plot_matplotlib\u001b[0;34m(obj, mesh, kwargs)\u001b[0m\n\u001b[1;32m    312\u001b[0m         \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"vmin\"\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvmin\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    313\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mvmax\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;34m\"vmax\"\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 314\u001b[0;31m         \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"vmax\"\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvmax\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    315\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    316\u001b[0m     \u001b[0;31m# Drop unsupported kwargs and inform user\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Fenics/comet-fenics/examples/beams_3D/plotting.py\u001b[0m in \u001b[0;36mmplot_function\u001b[0;34m(ax, f, **kwargs)\u001b[0m\n\u001b[1;32m    202\u001b[0m             \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcollections\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mLineCollection\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    203\u001b[0m             \u001b[0mnorm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mNormalize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mC\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mC\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 204\u001b[0;31m             \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmesh2segments\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmesh\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    205\u001b[0m             \u001b[0mlc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mLineCollection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmesh2segments\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmesh\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcmap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'viridis'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnorm\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnorm\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    206\u001b[0m             \u001b[0;31m# Set the values used for colormapping\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: Matplotlib plotting backend only supports 2D mesh for scalar functions."
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAADnCAYAAADl9EEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsvUlsI+l5Pv4U930ntVD7vu/qGRgxjAAe5B7YCYIsSGCkJ0ESxJd4Mr7kEsRwI0GQUzDjAEFu+dvjHJxDDmkEyQ/wkpZa3dqo1kZRlERx3/el6n9Qvm9IiSKrqJJFe/gAg2lJxaoiWe/3vcvzPi/DcRw66KCD9oPkqW+ggw46qI+OcXbQQZuiY5wddNCm6BhnBx20KTrG2UEHbQpZk793UrkddPD4YOr9srNzdtBBm6JjnB100KboGGcHHbQpOsbZQQdtio5xdtBBm6JjnB100KboGGcHHbQpOsbZQQdtio5xdtBBm6JjnB100KboGGcHHbQpOsbZQQdtio5xdtBBm6JZV0oHjwiO41Aul5HP5yGVSiGXyyGVSiGRdNbMDgCmicBXp2XskcBxHILBICqVCrRaLbLZLC4uLjAxMQGJRAKZTNYx1i8O6raMdXbOJ0C5XEa5XEYymQTLsiiXy0ilUlAqlWAYBhzHoVgsolgsAkDHWL+g6BjnzxHEjS2Xy+A4DqFQCAAwOjqKVCqFZDKJYDAIvV6PZDIJk8kElUpV11jlcjlkMlnHWH+J0THOnxNYlkWpVALLsuA4Dq9fv4Zer0dPTw/UajUcDgfi8Ti6u7uRz+dRKpVwdHQEs9mMrq4uRKNRmM1mKJVKcByHQqGAQqEAoGOsv6zoxJyPDI7jUKlUUCqVAACRSAQWiwWFQgESiQSVSgUymQyFQgFHR0dYWFi4c45isYjr62vEYjHodDoMDw8jFArBbDZDoVCA4zhUf48dY/2FQyfm/HmD4ziUSiVUKhWwLIvDw0MUi0WYTCZoNBr4fD7k83k4nU5qUPWgUCgwODiIwcFBAKAZXpfLBaVSienpaUQiERiNRsjlcnAch5OTE6jVathstppMsEwmA8PUfRY6aDN0jPORwLIsisUi3dGi0Sj0ej36+/upcZDkD3BjgNPT07zOLZPJMDQ0hKGhIbprptNpXF5egmEYLC0tIZfLQaFQQCKRgGVZ5PN5+npirGRn7Rhre6Lj1oqM6qQPAFxdXSGZTGJ2dvbOsfF4HPl8HmazGeVyGScnJ7wNtNH1GYaB2+1GLBYDx3FYWVlBOp2GWq2GTCa74wZ3jPXJ0XFrHxskq8qyLBiGwe7uLiQSyb0GZzAYoNPpUKlUAKBmd2sVxLD0ej2sViv0ej0kEgkSiQSOj4/BsiwWFxdRKpWgUCggk8k6O2ubomOcIqFSqSAYDCIajcJqtUKtVmNsbAwajebe14TDYSQSCQwMDIh+P5lMBhzHwWg0AgD6+/vR39+PSqUCiUSCUCgEv98PjuMwMzNDa6lSqbRjrG2CjnE+ENVuLMuyCIfDiEQiWFhYaGiY1a8HbuLIeq6v2JBKpQAAp9MJp9OJcrkMiUSCYDCIq6srcByH8fFxqNVqSCQSaqwnJydQqVQ1CaaOsT4uOsb5AFTXLgGgUCiAZVm8//77vMoXUqmUGgvHcbi+vqYZ2YfC4XDQczeCTHbzCHR3d6O7u5vGyvF4HB6PBwzDYGhoCJVKBRzH0QRTLpejRtkx1sdBJyHUAurVLr1eL5aWlsBxHH3g+Z6LxKlv377F6uqqKPeYTCYhk8l47d6NUCqVwHEc/H4/rq+vIZPJ0N/fD5vNRo2V4zgaZwMdY20BnYSQGKiuXVZnRefn51EoFBCJRNDf38/rXPF4HNFoFE6ns6asIgai0Sg0Gs2DjZPUXnt6euB0OikXOJfLweVyQSqVoqenB93d3TRT3NlZxUHHOAWgunaZy+XAsix6enowMjIChmGQSqUQj8d5GyfLspSCxzAMhoeHH/P2H4SLiwsYDAbYbDZqsOvr65TzWyqVsL29DZlMBrvdjr6+Prqb3jZWmUxG/+sY6/3oGCcPVLuxDMMgEAjA7XZjbm4Oer2eHtfK7ld9vJg7p8FggEKhEO18ZFe8DYVCQa+zvr6OQqGAfD4PlmXx+vVryGQyWK1WDAwM0PdXqVSQyWRwfX2NoaGhmo6bjrF+jo5xNkF1TAjcxGCRSATPnj27Q7dTqVSCEjpGoxFarZb+fHZ2BpvNJsp963Q6UR9yhULBK8GkVCqhVCoBfG6s6XQaAPD27VswDAOz2Qy73Y5sNkv5xSQRBYAaq0wmg0Qi+cIaa8c4G6DajU2n03j37h1WV1cxNzd372tIWxcfFAoFxONxWK1W0R9An88HjUYDh8Mhyvn4uuq3UW2sy8vLyOfzSCQS4DgOqVQKOzs7MJvN6O3tpewlYqxkt/6iGmvHOOuA1C6Pj4+h1WrBcRw8Hg/m5+cb7h6lUgkXFxew2+28rlMqlRCNRuluyfd1T4GzszOYzWaYTKYHnUelUkGlUqFSqWBpaQkMwyAWi4FhGBwfHyOXy8FkMqG7u5t23FTvrMRYyX+/zMbaMc5bqK5dsiyLdDqNvr4+XnXDh8acPT09Ld1zPSiVSkElnWYgO5mY58tms7DZbFCr1QCAsbEx5HI5xONxcByH8/Nz2nTucDhoL2u5XKZlrF9mY+0Y5//hdu0yHo/j8vISPT09UKlUvM4hk8kEuZFqtRpOp5P+vLOzg/X1dWE3fg/ENHTg/oRQqygUCgiHwzUxNsMwNeWfgYEB5HI5xGIxlMtlRKNRhMNhmM1m2Gw2qhLxy2qsHePE3U4ShmFwdXWF999/n67qfCCVSmG1WnkfTx6aUCgEnU4n7KabwOv1Qq1Wi+Yqj46OivqQsyzblEV121g1Gg0MBgNisRjy+TzN+JrNZlgsFoTDYbpA/jIY6xe+RZ4kfcrlMgqFAl6/fo18Po/5+Xmk02kEg0He56pUKtjd3eV9fCaTgdfrRSqVwtHREa15hsNhQYmleiByKGIhEAiI0jVDQHi6QsAwDLRaLfr6+mAymWCxWGhPayqVQi6Xg8/ng9/vR7FYpGWZUqmEXC5HdZoymQwKhQKlJLYrvrA75+3aZTKZxP7+PqampuhuWSgUalL8zSAk5gyHwzg4OIBWq8Xw8DB9LXnQvF4vbZxOpVLQaDSixpBCkUqlRN3dpVIp73DhPjAMA51OR+8rk8nAbreDZVlEIhE4HA6cnJzAZDLBbDZDrVZThtftnbVa2bBddtYvpHFW1y45jkMsFoNer8f6+npN4Z7wRvmCrOzNEAgEKBc3FovR329vb2NpaQnDw8MYHh6mDJtIJIKjoyMwDEPLESqVqqFbKBaBnkDsHSaZTCIej2N0dFS0c95mWHEch/7+fsRiMUp4ODo6gtFohNlspjHrbWO9rb/0VMb6hSO+V9cus9ksdnd30dvbW/dhJgb80BWegCRBuru7qeFXPxQbGxsNE0LlchkymQxnZ2cIhUKQy+W0zUwul9c8RIFAAEql8sGlD4JUKgWVSnWvzpFQhEIhpNNpUSmLl5eXMBqNNaytahCvJBaLIZvNYnp6Gqenp9BqtTXKhmTRZhiGGivDMMhms4JyCgJQ1/q/MDEnMYRqLuvFxQVmZ2fv3WVKpRIymYyg67x69aru76PRKDY3N6FUKmksVCgUsL+/z/vcxK0dHh7Gs2fPMD09Dblcjuvra7x69Qq7u7vI5XKoVCrIZrMPjlurQbpTxIJCoXgwKf82MpkMZXLVA8MwMBgMGBwcpOoUNpsN+Xweh4eH4DgOl5eXCIfDqFQqkEqlqFQqyOVyODo6wh/8wR+Ier/N8IVwazmOQyaTwfHxMcbHx3FwcACn09lUryebzSIWiwlaLW/HqNUCX6urqzW7cHWMynEclpeXeV8HAD3X4OAgBgYGkM1mqbFeXV3RTKdWq32wa+b3+zEwMCAaX9dgMMBgMIhyLgI+GeDbMBqNVC2C3FcsFoPP58Pc3BxOTk4A3OQI+IQsYuKXfuckMpIcxyEej2NjYwM2m42XwQmNOW+jWCxia2sL0WgUY2Njd9xjhmHoCAYAOD8/b/laJN6VyWTo6+vDysoKhoaGIJfL4ff78erVKxweHiKVSrV0frHrnH6/H5eXl6KdD7hZpB66GxuNRgwNDWF5eZkmiaLRKD7++GP813/9F/7sz/4MLper4Tk++OCDptdxu9148eIFXr58CYZhvsUwzJ3445d257xdu4zH4wCApaUl3l+gRqMRnO5///33AQC5XA5bW1sYHx+/dyFQKpWYm5uj7mc0GhUtQVIqlSivtaenB11dXUgkEpBIJIhEIjg7O4PJZEJXV9e9MVo1rFaraPEmIL6xAzdxMTEosaDVamG1WvE3f/M3+Pd//3f8+q//+r05iJcvX8LtduPly5dNz/v1r38dr1+/Jj9+CuB7AL5efcwvpXFWU/BKpRL29vag0WiwuLgoaGWVy+WUtM0X7969g81mg91uv+PG3kapVMLBwQGmpqYEXYMPYrEYdWmBGy/AbDYDuHngjEYj4vE4KpUKkskkTk5OKPOmnrGaTCbRSzliK9EHAgHREmAEhMGVzWZhMBjwq7/6q/ce+9WvfhUA8OGHHzY859bWFiwWC/2Z47g4wzBfvX3cL5VbS3ZLouVDmD4DAwOYnp6mg4P4IpvNwuPx8D6+VCrB5/MhEomA4zheWd7qwv7Y2Jig+3sIZDIZbDYbTCYTDAYD5ubmoNFokMvlkM/nsb29Da/XS9u9jo+PaTJNDPT29qK3t1e08wGokUoRC2dnZ0gmk8jlcqLVed1ud71FJMowzEr1L35pds5q+RCO43B6egqVSoWRkRF6jM/nw9DQEO9zki5+vri6uoJcLuctDH37QSJlFTFgMpkEuaEKhQJdXV3054mJCcRiMcRiMajVaiSTSQQCAXR1dQmiNN4Hv98PuVwuammip6dH9N2dLEjZbFaU9w3chC988EthnNW1y0qlgq2tLdhstpZ7EAnkcnnTeIzjOFxcXKBcLmNkZESQWyWVSrG0tER/Pj8/F63/kkhbPuT15GHkOA5ms5mWGsbGxnBycgK9Xk/rg0JRKBRE3+VUKpXo5yQllWw2K1rpx2Kx0BxI9a9v/+IX2q29XbuMxWJU/7UeUbt6F+UDpVLZkGnDsiy2t7eRSqXocdlsVtD9X1xcCLonviBTycQAwzCYmZnB8PAwxsfHAdxIb+bzeRwdHYHjOHi9XoRCId50x8dICB0eHop6PgCYnJyEWq0W1a0dGRmpu3tyHLdV/fMvrHESCh4Rc3a5XPB6vahUKvd+iEKzjYVCAW/evKn7NxKLDQ0NYXZ2lmYIhcSoAGqI9dVuZbthb2+PZpUZhqElh/n5eTAMQwf+bm9vo1KpIBAIIBqN0lETt9Hb21uTFBEDhNEjJq6ursCy7IN3TrfbTXfLlZWa0BIMw4wAuJPi/YU0TjKvkiQAyDTopaWlhjGH0JW1Xp2TuHa7u7t0nF+ruE2UF0s/CICoVDvgLrniNsxmM0ZHR7G6ukqlMMPhMLa2tpDP5xGLxZBMJun7TaVSojKYAIieqQVuNIkB8DLOra0tvHjxAgDw0Ucf1ZRUPvroI3z/+9+nP3/ve9+jdU4AXwPwh7fP9wsVc1ZnY4+Pj+F0OpHNZjE5OSnoHHxX13orcTQaRSwWw7Nnz+rW0/r6+njfC3CTeCHY398Xrdm6u7tblPO0CovFUrMzZjIZ+P1+pFIpTE9P07KHWq0WbbcTk0RPQGLOXC7X1DhXVlawsrKCb33rW3f+9oMf/KDusQDAcdyLeuf7hTHO6tqlVCpFOp1GIpEQNDLPbrcLMk65XE4V2NPpNC4uLjA1NdUww2ixWHhfg2GYO26fWLGYx+MRVeBrYWHhQfdltVphtVpr5DH9fj8uLi4wPj4OjUYDjuMelBF9/fq1aIsbAcld5HK5Dn3vNm7XLlOpFDiOQ09PD+bm5gSxQUjfpJBr7+7uwufzYWdnB319fU1fv7e3J6j8cnp6Sv9NxMTaEWdnZ4Le130g3kh/fz9mZ2fx7NkzmEwmSi5/9eoVQqEQSqWS6G5vKyDufCaTEZ2o3wxtvXNW1y6Bm90gFAphcXGRZg2F4O3bt5ifn+ed+ifsme7ubjx79oxXDa0VkS/ympmZGUGva3ZOMUHkLMVCda8k0bI1m81UWC2bzeLk5AQsy1KBNZZlG34HYtL2CDweDxwOx5PsnG1rnNW1S4ZhEI/HUSqVsL6+DolEgp/85Cf40pe+JOicZEIWH2QyGezs7KBcLgvS4RHaadHV1YVAIACZTIbz83MsLCyI8pA9RrO1mAZ/eXmJvr6+O0kriUQCiUQCg8GAlZUVKouZy+VwcHAAiURSM5ulupZ7OwsqBsh77hgnagnrDMMgGo3C5/NhYWGBckNbBWmmbYZisYjt7W3Mzs4KZpwIHRvPcRyi0SjVugkEAiiXy7BYLA9q9fL7/VAqlQ/+zAjGx8dF5cLyNXYyokGpVGJtbY3OZSHfkUKhgMPhQFdXF46OjgR//s1AyluZTOaLHXNW1y4ZhsHp6Sk8Hk/dbGwrrJTZ2dmGcQMZEiuRSPD+++/DaDQK5uNubW3x4qCm02l4vV4EAgFMT09jYmICMpkMRqMRMpkMHo8H8XgcsVgMfr9fMK81n8+LSgcUQq7gA71e31KpR6FQQKvVQqlU4tmzZ5iamoJGo0G5XEYwGKQ5AjHAcRz1mgqFgmiKGHzRNsZZqVRo0oeQr7u7u7G6ulrXEFuZY3l8fHxvP2M2m8WrV69qtGMAtPRFN3Od/X4/HUNQvXuQdrbe3l7Mzc1Ralw+n4fL5UIwGEQmk0EkErm3uP9YEJvJ1NvbK0rjtkqlovVNo9GIkZEROtZha2sLLpcL19fXLcXL5XIZ7969oz+L3UXTDE9unCTpQzJzgUAA29vbKJVKDYfxvH37VvADWiqV6hbTiazH1NQUhoaGHhRbNWrQJto0+Xwe6+vr0Ov1NTu5x+O58540Gg1t/iVlkWg0iq2tLVxdXaFYLNYU9wl6enpEc2kfA/v7+6J2ucjlckxOTkKr1cLhcFAxtL6+PjpLdX9/H4eHhwiFQrxyD6Rs91QZ9CeNOVmWhd/vh1QqhdFoRKFQQCgUwvr6elOXh2RxhSRPbieEWJbF0dERANzbUylUgKqayF6NQqGA7e1tjI+P13TGLC0tUfczmUw2fWi0Wi3NVBORssvLS6RSKdqGVSqVkM/nqeSjGOjt7RU1ISR2gqlUKiGRSNTUSYlmEEnSTU1NIZlMUukZIkFisVhgMpnu7Iwcx9HF8zG4wM3wpMZZLpeRSqVQKpVwenqK5eVlLCws8HqtkMwrwejoKP0COI7DmzdvYDabGxqg0NjW7Xaju7u7ZkckqgiTk5N3drONjY17DboZiDTJzMwM7cgpFos4PDxEOp2GxWLB9PQ0KpXKg1upqnV2xEB16CAGisUiYrFYQ2aUVCqlJRvghhcdj8epPhCpr1osFuj1eqhUKoyNjX0xd07C9AmHw1hbWxP0Zc3Pzwt+4OLxOBQKBYrFIrRaLRYWFpruLO/evRNUskmn09Q15TiOzgO5TxWhOmkzOjrashFVjx1YWlrC2dkZ5HI5SqUSdnZ2wDAMurq60NfX19IucHBwICr7ptEYxVbQirgXaTgnnGaHw0HFvfr6+hCNRhGJRCCXy3/u8SbQBjFnV1cXlpaWBLfjBINBwTFLMpnE6ekpzs/PKTmbD4SsnGRHJ6MZAoGAIFUEMVg4wI2rZjaboVAosLa2hoWFBRgMBrAsi42NDezs7IiW1WwFLpdL1KSWSqV6cFcPaTifnp6GXq+n4mt//dd/DY/Hg9/6rd/C//zP/4h0x83x5HVOg8HQEk0rkUhAr9fz5mKS+inDMFhdXeW9cwjl446MjEChUODq6goWi6UpEf7Zs2fU+K+vr6HX60Xp5r89iVoul9PFaH19HblcDul0mrr3arUaFouFJlNug48ImBBks1nRYzgxu3AA0NLWP/zDP+CP//iP8fHHHzfMcbjdbnz22WdYWVnB1tYWnj9/fm+nzNbWFjY3NwEAH3744bcAfMZxnLvm+uK9ldaQTqcRCoUEi1yRceV8EIlE4Ha7MT09LbjnTygf9+rqCuVymXcx/PT0lDaBt0L9uw+BQOBe4vvt6V1kHksikQDDMHC5XFAqlbBYLDAajZBIJKKLkImdYEkmk0ilUoIb6hvBarXCYrHA4/FAp9M1zYdUK+qtra3hD//wD+90oxC8fPmSdq88f/78BcMwnwCoUQZ7UreWYRiqqi0UTqeTF2PD6/Xi9PQU8/PzLenQvnnzhlcxn+M4HB8f0x2TLyKRCH1IHQ7HkwwrkkgkMBqNGBgYAHDTxmYwGCgB4+zsDD/+8Y/pTisGSJO2WLhN5RMDREOJT7vYbUU9k8nUUCLzk08+qSdVUoMnjzmVSmVLyt/12q2qQep/NpsNa2trUKlUSCQSghk/fLLCZKYKGZcu1MDIQ6rX60V7wDQaTctFfplMBrvdjomJCcplBW7qsIVCAX6/H9fX1w+qUwYCgZZfWw8KhUL0GafZbBaFQoGXuFc9RT2LxYKtra26x3/00UcYHh7Gp59+CoZhngP46PYxT+7WqlSqlhqDw+EwlEpl3RUtFovB5XJhcnKyxvBbKb80e8ATiQT29vYwPz8Pu90OiUQiqCdxamqK7kZnZ2cYGhoShcMp1tBc4OY7UigUNMOq0+kQDofhcrkwNjZGpWJMJhPvunMwGBSVnK/X6x9F3Esul/NSQeCrqEfw/PlzxONxfPLJJ8CNO/sSQM1W+uRuLSEwC8V9hsayLM7Pz7GysnJH9kOj0QhObMzPz9+bac1ms3C5XFhaWqKLgFqtFrRzismSqcb5+blgL6ERqssoOp2Ospb0ej0dWUBGT2QymbqspcfE9fU1/H6/qOd0OBwwm83IZrNNF8x6inqNDPbFixf41re+RWLUTwD85+1jntytbTXmNBgMNR9YsVjE27dvUSwWsbS0VHf30mg0gnVmjo6O7pC+WZbF2dkZVCoV3nvvvZr78Pl8glTvzs7O6EOs1WqfpJ7GB430lwwGA8bHx7G+vg6z2YxyuYzLy0vaOF0sFpHL5WpeI0Q/mA9aqXM2Aylt8THO+xT16rWxvXz5sub3HMd9CuCz26LST/4kCKk3VkOv19MPLJlMYmNjAz09PQ3riWTsgBDcnm6dy+Xw6tUr2ih8+4EQmnGtnmg9NDQkmnCx2ENfk8kkr+OIMt/MzAyePXsGq9VKJTRfvXqFQCDwKOPeJRKJ6Mm0UCiEQqGAXC7X9Hu5bYRut5uOZyA/k531vlj0tjTmk8ecDMO01CQbiUSQTCbR09MDuVyO5eXlpnFBKzFndYaXTMGuR8Mj0Gq1ghIxZBYHcMNGcjqdotQUSeb1KUHKVgaDAYuLi5Sckcvl4HK5cHFxAafTia6urgdnW4UKq/FBpVKBRCLhtXMCnyvqkTrn9773Pfq3jz76CB988AGeP3+OlZUVuN1ufPrppwBonfP/u32+J59sXSgUsLW1JdhA/X4/jo+P4XA4eKvvZTIZXF1d1SjeNUOxWIREIoHX60W5XG76WvJ58t21MpkMFAoFKpUKjo6O0NvbK8rcyuvra1qrFAOpVEpUIgLhFJdKJTAMg93dXcryIcYqZOf3er13Zm0+FMFgECaTCf/4j/+IgYEB/N7v/Z5o576F9p1sXT3Mhy+urq5gMBgEyWJqtVpBhgncBPUkluUzaMjr9eLq6or3+ff29mjMLaYbetsdfyj4urV8YbPZIJfLodFooFarsb6+jrGxMSiVSpRKJWxsbODg4IB3UovMYBUTZN6pmHNShODJjVOoGp7X60U4HMbKygoWFxcFXSufz2Nvb4/38el0GrFYDF1dXZiamuLldgmNOYnCn8fjwdjYmKgsITEhJg+X47gadx74nLVEBjCtra2ht7cXlUoFLMtia2sLp6en9ybbHiMhdHR0RHt9f94SJUAbGCfAT9WgXC5jZ2cHyWQSZrMZmUwGBwcHgq7DMAzvXfr6+ho7OzvgOE5Q25iQRAzHccjlcujt7cXAwAAuLi7gdrvx6tUrmsVtVR7yMcYdiAWO45qWzwhrqbu7GxKJhBL3E4kEANDxG4S1NDo6KjoJgRj8U4h7AW2QEAJuXNTBwcF7C9gk7uvq6qKEBYZhBNcI+SaE0uk0AoEAnj17hmw2KyibzDcx4ff7YTQa8ZWvfAUsy6JcLqNYLGJwcBBGo5GqNhAdXIfDgYGBAd6xWDqdFpU1I2bCpZXkD2EtEXIFKV0Qoe+DgwOYzWY4HI6W9KXqgZS2xJwwJgRPvnMShb374qOrqyu8fv0aDMPUMIlaqY/KZLKGiadCoYC9vT1otVo6d0UoHzcUCjV0ATmOw9HREXw+H+RyOba3t2sGBBHjUygUkMvlWFlZwdLSEtV03dzc5NXulUwmRRXlEnvneKhYlkqlQm9vL21mIBpULpcL5XIZgUDgwVpLU1NTYBjmi2ucwP2G5vV6EYlEsL6+fmdXVSqVgmNO4KbeVA+xWAybm5vo6emp2ZnC4bAgUgFRp78PmUwGEokEy8vLkMlkVFcIuNmd6mVEZTIZ5d2ura1hbGyMili9efMGh4eHCIfDjxqrijleTyqVYnZ2VrTzATde0cDAAP1clUolYrEYtra2kMvlEI/HBbOWSEuXmOP/hKAt3NqhoaGa2iCZg9LX14f+/v66bhzHcTg7OxOk/M4wDMLhcE3GlvR5siyLlZWVO1k5oTvnfdOw0+k0Dg8Psby8XJP1rX5vhEXTqJh+u91rcXER8XgciUQCVqsVh4eHUKlUUKlUou12Yht9sVjE6empqBqz/f39NZ+byWSqYYNls1n4fD46SEkmk4FhmIZZWPK++XSlPAbawjiJWyKTyXB9fY2zszPMzc01jEsYhkEoFGppLANBuVzG/v4+9Hr9vX2AOp1O0MNptVrvEBRCoRCOj49p21o1lpeXqap5LBYTPMxHIpHUTPQaGhpCNBpFMBiEwWDA5eUlZDIZLBbLg6QoxawfkqFUYqJZLH57kFI0GsXl5SUKhQLGx8dpfFlt4OR8TyEoDbSBcTLoFT9XAAAgAElEQVQMg+vra6rlks1mec0labUmSOqiRNfU6XTeSetXw2g0CjJOQvciCgr5fB4qlQpra2t1jePs7IxmJMWAUqlET08PCoUC0uk0jEYjwuEw9vb2MDs7i3Q6DalUCoPBIOiafGq8fPEYSnZut5tXdppclxgrCSvi8TjOzs7AMAyGhoZgMplofqJQKIiWZBKCJzdO4ObLevfuHZaWlgTNWGylO79UKtFhuysrK00XgUAggEKhwFsiM5fLIRaLwWw2Y3d3lwod3wfSc6pQKGCxWER/CPR6PfR6Pb1/Mifz6OgIs7Oz4DgOUqm0qXu3tbWFtbU1Ue5JqVSKrqwAtLZgkwWKGGupVKLG6nK58KMf/QgMw+D4+BiTk5M/V3nMJ08IJZNJhEIhjI6OCqaHCS2lcBxHOZ1Eqa4ZhPJxSYy6vb2Nrq6uprIZ1V+2SqUSjbxNRhbchsViwdTUFNbX16HRaJDL5XB4eIiNjQ0kk0kUi8VHmxlKUCwWEQwGRTsfIF7/qlwuh0KhgF6vh06nw+///u8jn8/j29/+dt2huI+JJ985dTodLxHpejg7O6Nd+s3AcRwymQwl2gsZoCu0P5MIP/NxG0dHR6FQKKjAtk6nE+VBs1gsDd8j+RupHRKDjEQiOD8/h0wmw8jICPR6Pd1dxUKpVEIymWwYTggF3+eAL0gOZGZmBkqlEv/2b//W8Hgh4l4A8Nlnn9F/f/3rX/8ax3Gf3T7myY2TNOqyLPtonRSJRALHx8dYXV0VLNhF3J1m4DgOHo8HwWBQ0ESuXC4HqVRK70mszOjFxYWgydbE+BwOBxwOB/L5PCQSCSKRCM7OzqDX65FMJkUh5T9GzLmzsyOqrq5CoaAhFp97FSLu9eLFC4yMjOBrX/saaSP7GED7GSdw4wq2kr3r7e1t+kX7fD6cn59jcXERDMPA4XAIejhisRii0WjDWJjQ7IrFIsbHx3F9fc2bOuf3+9Hf3w+VSvUg3R+xQUgCdrsdJpMJ7969g0QiQTQaxdnZGSwWC+x2e0v1P7Vajd7eXrFvWVSUSiXKDmvmNQgR94rH4/jOd75Da+cmkwkcx9Xlrz55zAncJAhaibUauX+VSgWZTAYGgwHPnj2jdard3V1BrJFm/NZcLofNzU1IpVJMTk4KHnxTTXR3OByilSzEHndQLpeh0+lgsViwsLAAjUaDQqGAVCqFvb09+Hw+3jkAsYSzqyH2eL5cLodoNMqLHSRE3GtzcxMjIyP47LPP8PLlS7x48QIMw9RNTDy5cTIMA6vViv7+fsGvdblcdXfcbDaLjY0NxONx6HS6mpVPKpUKTvDcdzxhoIyPj9PFRS6XC3L9HA4H3S0vLi5EU6Xr7++/o6HUKm57GnK5HF1dXbBarVRPqFwuU+WA4+NjRKPRexfBTCYjekJofn5e1PORCWN82sWEiHu53W5sbW3hq1/9Kr761a/i+fPnQB39IKANjBO4iQlboYfVMxyO47C/v4+pqam6CQehbli9nlEy3UupVGJ1dbVm1VSr1YISHTqdTnSlcuDGnReqCHcfZDLZvX2zDMNAp9NhYGAAfX19kMlkMJvNCIfDcLvdYFkWFxcXyGQyNYoSYsacLMtiZ2dHtPMBN9+jzWbj1ZEiRNxrZGQEIyMj9Jn5v/+P1Ns928I4iQqfUFTvUBzHwe12I5PJYG1t7d5M2dTUlKC4rlgs1pDMK5UK9vb2cH5+Do1Gc8edymaz2N3d5X1+j8dDB/qKPXVLrGZrlmWbCiATSKVS2Gw2TExMYHx8nHaguN1u7O/vAxCfq0q6esQEGXWfyWSa7pxCxL3uKa3V/XDbwjhbVeAbHx+HSqVCqVSiyuwajabhqry/v490Os37GpVKpeaDd7lcMBqN9xbRWxX4Am6I761o+D42iMvaCqRSKZxOJ+bn5zE7OwuWZZHP5+H1eulul81mH5SlfoxG63A4jEgkwotXK0Tci+ya5Of/+7/79pwUoA2ytYTI3UrMcHBwgK6uLto+xOfBJnMshdwfEfZiWRazs7MNHwShKnBECgO40ayRSCSixIpOp1M011EsN5QIfun1esjlctoj6vV6kUgkoNPpMDMzQ+M9vpDL5Y8yUpBIlIgp7gUAP/jBD/Cd73wHo6OjOD09BYCv1zvnkwt8ETUAt9stmMS+ubmJQqGAL33pS7wfHpfLJWgke6VSwcnJCWKxGBYXFx9FS+bq6gpGoxFerxcKhUKUem8kEqFMl4eiXC4jHo+LlmDy+/0oFos175PE8VqtFi6XC5lMhg42bqYuQSaii9kQHgwGoVAo8OMf/xjb29v4zne+I9q566Dum3vynRP4vJVLiHGenZ0hl8thfHxc0KpOGmj5gOM4Kipcr6e0HkqlEuWtNgPLsvjZz34GqVQKk8mEfD6PYDCIcDiMubk5OhC3lV0rlUq1pHBfD6T5WyxIpdI75yNTugFgZmaGLggSiQSHh4coFouwWCx1Z9EQxpGYsNlsYBjmydrFgDYxTiE9k/l8HtlsFk6nE4ODg4IfXL/fD61W27SeWCgUsLOzg+npaZRKJUFuViaT4XVcPB6HVCqlEi3VfaZyuRzn5+cIBAIwGAx07ufPk3hNQGJEsVxHPvREMnUauFlQs9ksjf3dbjcdD2+xWB4l5jw9PaXZ2i+scZKHjU/jbTQaxcHBAaanp6FQKBAIBCCTyXjR6wiy2SykUmlD40ylUtjZ2cHk5CTUarWguiifhFAmk6FjCa1WK/0MiLwjubehoSEMDg4imUxCKpXC4/EgFovBarWiu7u7YQeLzWYTlQ8r5qLg8/kgkUgEJb+qG8yHhoaQSCQQjUah1WqRTCYhkUiQSqWg0+lEuddKpULrnE+hggC0gXES5HK5hnEgx3G4urrC6uoqLV8UCgXBtL9mXSakfkkU5DmOEzTOXCqVNpRPCYfDODw8pPMpx8bGwLIsisUiMpkMyuVyzcJBxhsAN4N8nU4nIpEIWJbF5eUlUqkUbDYbzGbznUZhsQyq2WImFJVK5UH3JpFIYDab6fNCxNG8Xi+lBRYKhQc1mJPZsdls9sHj7FtFW5RSGIaBx+Op+7dyuYzt7W1kMpk7E7+ETLcm6Onpqct7ZVkWBwcHOD4+pmLH5N6Eljeur6/r/p7U4tbW1miNVihZQKFQoKenh/JTe3p6kEwmkU6naUdJOp1GMBgULQ4jDdxi4TFa0ABgdnYWZrMZSqWSahSHQiFks1mabeeLyclJaLVaXnNSHgttYZwEt93BbDaLV69e3Uuwttvtgle1UqlUlwN6fHwMhUKBhYWFOw/OmzdvBF3jtnGyLIv9/X14PB7YbLYadzSXy9GHy2AwCJqCJpFIYDKZMDo6CpPJBL1eD4VCAY/Hg0wmg3w+j1Ao9OACfTqdxtHR0YPOUY2+vj7e3TJ8wLJszXem0WgwNDSElZUV2O12sCyLQCCAzc1NXF9fo1Qq3Zl6dhtHR0eU/P6Fd2v7+vpqVtRsNguZTIb5+fl7M45E6U4I6TmZTKJSqdCdK5VKIZPJCGrzaoTbhs1xHN6+fQuz2Vx37B05nmRmHyLlSHbVnp4ehMNhsCyLZDKJ8/Nz9PT0wG638yJqPDZCoRDUarUo7WdA85moOp2ODikmJZvT01MUCgX09/fT2SzVMXo6nX5SWUygTYyTkN+Bm1Xw+PgY2WwWS0tLDWMGksETomxe3Z7m9/vhdruxsLDQ0DCFijuR9rJ0Og2ZTIbZ2dl7kzfVrVPxeBy5XK6uEQuFyWSiLXLA53VEt9uNXC5H2Uh8G6nF5P9mMhlRk1VEbbAZSByu0+mwuLhIhdXI1DOZTIbe3l44HA6Uy2Wq9v6FNk7gRhd1enqaSogsLS01Xd2FdpgAN5o6SqUSHMchkUjwUmEQqo8rlUoRCoVwdHSExcXFhm5RLpejEo1i7maXl5c1zdakjjg/P0+5qPUeSnJsNYgOkVi47YY+FKTDRaiyglQqpRza9fV1FAoFFItFFAoFsCyLf/mXf0E4HBZdKZAv2ibmrFQq8Pv9GB8fx+joKK8vT6PRCGatSKVSnJ6eIp1OY3JykteOsLGxIega7969w/n5OdbX15vGK4lEghLflUrlz2WVlkgkUCgU0Gq1WF9fx/T0NO3P3NjYwOHhYU2SKpVK3ZuwawU9PT2iG7sYIYlSqYRer4dKpUJPTw/m5uYQiUTwzW9+E7/xG79x7+vcbjdevHhB+zP5Ngl8+OGHDf/eFjsnaSkaGRkRtKLK5XJBmTTSGK3T6QQ9HHxXTuKSMwyD1dVVXu+lWp5EzB1KJpPxdh2r3cK1tTUkEgna1bK3tweNRiPq7lEsFkWdQk0U3sVEMBikImg//OEPG8bHQiRKCF6+fIlPP/0Un3zyyb3HtMXOqdPp8OUvf1kQmQC4KdrzHSNPhtSOjo4+Smq8WCxic3MTSqUSdrud9yJjNpthMBjozBixdqi+vj7BnyfweQ2xq6sLMpmM0h1zuRw4jsPJyQnC4fCDElfBYLClmaz34b7ymBgg/Zz35T6ESJQQxONxWCyWppn5tjBOi8UCt9stuC2JT6sZeaAODg5oQV+ocT579qzpNTKZDIaHhzE0NCQooaNQKB5FN+jy8hKRSOTB51GpVBgZGcHCwgKAzxuLSW+m3+8XPDBJbLrd2dmZaI3lBIODgwBuFt1G348QiRKCly9f8prk3hbGCbRGKJDJZE1jOo/Hg3K5jJWVFUgkEiqrIQSnp6f3UvKCwSAtlRDOKBmAwwehUIgmNMTU/SmXyw/a3aqRz+eRSCTAMAwsFgvGxsawsLBA2+9OTk6wubkJjuOQSqWaXtdqtYq6IJXLZdGVFUguotm4QqGLwsuXL2t6PRuhLWJOoLXMq0wmu3eMPCmzDA0N1XxxmUwGJycngjKwhF1yO4bzer0IBAItTTsjqObiVlPS2gn5fB7JZPIOYZ1hGDrOgtSoQ6EQbVebm5sDx3F34kvSzykWxN6JK5UKPB4PL0KIEIkSMjKCL9GkLYyTYZiWOkw4jsNPf/pTfOlLX6r5PeGvkrar29cSugjc7pphWZZyN/v6+u48GEJIEWq1mjJ4kskk4vG4KP2cfX19P1eiAbkW0cjJ5/OQSqU4OjpCPB6nJAyZTIbT01OMjIyIlpkW2jbYDNXG3oyjLESiZGtrC9FolHpW8Xgcn376KT788MORtlRCIIhGoyiVSoJ4rPXYOMDNm15fX6/rOgkd6QcAExMT9MsqFArY3t5uOACJTzxBQFqeiBsqdMTEfYjH41AoFKKwcLRareCdjixQExMTqFQqtD3u+PgYsVgMgUAAfX19ouygl5eXcDgcosljEvedD/hIlJDd8mtf+1rNsR9++CGeP3+O58+f1x0a2zYxZ7lcFqTtQ0CMplKpYGdnB5FIBGNjY/fGNCqVSpDxADdJATKNant7GyMjIw0L3ltbW7wXAL/fj/Pzc0H3wwfpdFq0jCjDMA9yG6VSKaxWKyQSCSYmJmgjM1HJJ1KarerZJhIJUWeIymQy9PT0oFQq8Sr5EImSly9f4rPPPrsjUfL973+/5vh4PI4XL14AQEPd2ieXKQFu3Ijr62tEIpF7JRgboVgs4vXr1+jv728qVcGyLGUj8cWbN29oW5ZGo2n6oP70pz/Fe++9x+uBvr6+RjKZhNlsRjAYhMlkgs1me3AM5fF4BI1jaIRQKIR0Os170poQkF01Eomgt7eXNgJYLBbeWXXSeytWrTORSMDv96Orqwu//du/jf/3//6fKOdtgLp+c9vsnEajsSVh6a2tLchkMszNzfHSkGEYhjeDA/i8THJ1dQWlUsnLaPgq8HEch+vra8RiMWg0Guh0OlxdXeHq6gr5fB7X19ctSYYCNx07YhHLxW7x2t7epnE22VUnJiag0+mg0+noAhqNRpFKpZruqsPDw6ImmKobrZ+KVwu0UczJsizVveEDjuPg9Xpp4zFfZo1Q6cpCoQC1Wo2ZmRneDwAp2zSDy+WiszMLhQLkcjlUKhWcTiftuNnb24PJZMLg4CBtX+JjKGJmMFUqleijHe57D2q1Gv39/XShTqVSCIfDODk5wfDwMAwGAziOq4kvxZ48TTSTOsb5fyiVSri+vubdn3l9fY1UKgWtViu4w4GPDGehUMD+/j4WFhYEGSZwo6bX399/b7ySz+fp0KNisYhXr15Br9cjFothbGyMKgX09vbSh7RcLsPr9SKdTqO7uxsDAwMNJSQjkQg0Go0obCitVivqQypkJyaURtLulUqlcHZ2hlKphKGhIVgsFlxcXIjaH0pKWjs7O0/WaA20iXEK6WXM5/Pw+/0YHBxET08PbzGtasRisYY7bTKZxO7uLqampiCTyXB4eAin08m7PhUOh+9NGJFzT09PQyKRQKlUYmxsDKVSCQaDAVarFXt7e7BarbTtC7hJfJF4vFwuo1QqYXt7m84tEXs+ZTVCodAdKcuHYHW17lCthiAlDaPRiKWlJVQqFVQqFaTTaWQyGezt7cHpdLZEWbyNYDAIlmV5jWJ4TLRNzCmXy5uOhYvFYnj9+jXlopL4USgT5vLy8t6/lUolMAyDpaUl+kULdYXv0yniOA7Hx8dYXFyk7lkgEMDp6SlMJhP6+voglUoxPj6OSqVC41yfz0dnjZD+SzIEeHx8HEqlEpVKBZubmzg5OUE8Hoder3/SVb8RSHPAQ0DkNY1GI8bGxugQ4lQqhY2NDZycnLQs00Ky85lMpuPWAjcPdKMPguM4xONxrKys1Dx0wWAQVqv1wV0O1bNWCI+UQKvVCoq56iUovF4vzGYzVlZWaDKEjF9fWVmh90/6LsmKzbIsNBoNTQ4tLS0hEonQxAkZkw7c9J0mEgkkk0l0d3fj5OQEFosFVqv1QQkTqVQqasIlkUiIdi7ghmxevcMtLy8jFouhUCggl8vh5OSEDkHmk9GtVnt/SuNsm50TACVTV4NlWbhcLvh8PgwPD9/ZDVrh5BJSczWIcHE9bdbBwUFB6nNkhyP/Pjg4QDweh0qlQrlcpmR8mUyGxcXFhgsLkZCcm5vD8vIyGIZBIpHAzs4ODg4OIJfLwbIs3amNRiN6e3txeXkJrVaLfD5PJ7j5fD6k02nBNUG73f6obvNDcXBwUPOzTCaD3W6H3W6nxP1yuYzr62tUKhWcnp4iHo/fmwHu7e2F3W5/UhUEoE12zvsoUuVyGVtbW7Db7fe6vH19fYJJ1NWrbD6fR6FQoKLN9XB6egqDwcBLDBm46U9VqVSQSqUoFApQKpV04la5XMbu7i7sdrvg3Z58RmNjYxgbG6N9ke/evUM6nYbFYqELTyqVQnd3N8xmM+3ekUgk8Hg8yOfzWF1dRTqdhkajaZpQCwQC4DhOtCFLrdSyW0U9T8RgMMDv98Pn82F6ehqBQICq9gGgerhPHXO2hXESVJO+0+k0tFotpqenGyZvWqFsHRwc4L333qNTmYlIdSMI2Z0lEgny+TxcLheWlpYwODhIM7CZTAYDAwO8Db0RyD3Pzc3RYr5CocDOzg6KxSLNipIdgoyKBz4fgREOhylJHUBdQxV7vF46nRatBgtAUIO6RCKhuyrwOU96f38fMpkMCwsLCAaDcDqdyGQyour1CkXbGCfDMHRF9fl8OD8/x8rKStMP/vLyEna7XVCWTiKRoFwu4/T0lIpHN7s3Ia6gRqOBy+XC7OwsNY5kMgmv11tXelMMSKVS6HQ6vH79GvPz8ygWiwiFQlAoFHC73ZBKpTCbzXQxY1kWfX19GBwcRLFYhFQqxbt375BKpegOLKZaQTWurq6aJv+E4L5xjHwgkUgwODiIwcFBStFMp9P413/9V/zwhz/E0tISAoHAkwhLt41xAsDPfvYzTExMIBAIYH19ndfDIbTVjBhZNpvlzbEdHh7mbVDFYhE2mw0Oh4MqxgeDQZyfn2NxcfHROkWy2Sx2dnYwMTEBpVJJ9XCAmxH0RHB6cnIS8XgcDMPAYDBQNkypVKJZ4kQiQQcIlctlmlRqR7Asizdv3rRUnrkNkvSbnp7G8vIyTRD+7u/+Ln70ox/V9dLcbjc+++wzOvrv+fPn95bctra2sLm5iXg8jo2NDXz3u9+9b5gugDYyTqJ6Zjaba+aHNIOQ3sBKpYLd3V0YDAZBblUoFALDMA1dUY7jcH5+jlgsBoVCAavVStvB1Go1VldXH20nIgTt2dnZup6GUqlEb28v3a1UKhUCgQDOzs7oQFuJRAK5XA6GYWAymcCyLMbGxpDNZqkiYj6fB8uysNlsD9K+FXPXbJUs3wiVSgUymQxyuRy/8zu/g1/7tV+791i++kHxeBybm5t0RufLly/xwQcfkPmcddE22dq9vT0olUrBsokOh4OXobEsi1KpBIfDgXw+L6hjg6TkG8HtdiOVSmFubg52ux1erxebm5vweDy8tWFbQTAYxJs3byCTyXjHXiaTCZOTk3j//feh0+mQz+dxdHRUE6sCN5+ZSqWCRqNBsViE3W6nbnImk0E0GkUkEhGcLRcz3nyMz/b8/JySEBqFPEL0g9xuN7773e/Sn9fW1momXtdD2+ycKysrLbX9XF5egmGYhqT5eDyOd+/e4dmzZ+jt7UUoFBK04lYLUd9GuVxGKpVCf38/7RW1WCy4urqCzWaD1WqFx+PBwcEBBgYGaMZTjAcqEAjg4uICy8vLLXNfGYZBV1cXurq6UCqVIJFI4Pf74fV6YTAY4HQ6EQqFMD4+DqlUSl12EksTEgWRGSU6sI3w7t07rK+vt3S/tyGXy7G0tCTKuQiIJ9HMOBvpB90OmVZWVvCf//mf9OfNzU2YTKaGrLO2MU6GYXB0dASn0ymotkSSO/chGAzi9PQUi4uL9AGWy+WCFgKbzVb3+Fwuh7dv31JCNsuydEeempqqyaaSMko2m4XL5YJCoUBfX19LWVuStDCbzaKO+iPhgdPpRG9vL/x+PxWcJnNeCKGBZVmo1WoMDw9TneFoNAqv14tKpYKxsTHo9fo7U6nF7LsEbrwaUgMXC6T3tFkpRah+UHV8+cknn9T0fdZD2xgncFNzFJq2l8vldd0q8gAbDIY7yaWZmRlB1yiVSiiVSncSAoeHh5iamoJer6cZ2YODA8zOzt5ZYBiGgVwuh1wux3vvvYdsNotyuYxMJoP9/X1YrVY4HI6mrikhZSgUinv1k8RAPB7H+fk5FhYWaOknHA5TpXgyBoMklYCb+H9hYYHWVCORCM7OzqDX69Hf308J7GKWJ8iCJya6u7tp2auRcQrRD6rGp59+it/8zd+8o4xwG21jnGQeotD4pV6KmxT69Xo9xsbG7vzd4/HAbDbzfkiy2SwSiQStw/r9fsjlciwsLND0O8uyOD09xdLSEi9Oa7XxLi8vIxKJIJVKQaFQ4PT0FHa7HRaL5c6ueHJyAr1eX5flJBaKxSKOj4+xtLREFyTSygV8XvMNhUIolUrQ6/VIJpOUqE9cXpLlzWQykEgkCIfD8Hq9dGK0GNzfx5hqvbu7i/X19aZurRD9IIKXL19iZGSElwJf2ySEgJs3K1TxPB6P3xFiPjg4QFdXV13DBIQP3SV1TsK/vby8hF6vR6VSoUNsS6XSHd4vX8jlcnR3d6O3txdyuRw9PT2IxWI4OTkBy7K4urpCKpVCKBTC2NjYoxpmMBiERCLB2travQQPqVQKu92OmZkZes+FQgEulwvBYBAMw6BQKNBFS61WQ6FQUFlNn89HF7yjo6MHSZSQcX+PgVwu11B6lY9+UPXOShJI5JjPPvus4fXbZucEPldlF2KgZHoWcNO1Ui6X66ruVaPZdOvb0Gq1YBiGtikRKUyO4yhvVSy1u9tTmyuVCjKZDA4PD2E0GmG322m3hNg10/Pzc0SjURpz8QXhr46MjNDv4/z8HNlsFn19fejq6gLLspQlpVQqYTKZwHEcrFYrwuEwMpkM7HY7vT5fyZFisYhsNitqBw45VzNBaeBz/SBS57ytH/TBBx/g+fPncLvdd2qxIyMjDV3bttAQAm5c0ZOTE1qT44tEIoGLiwuYzWZ4vV5ebmWhUIBUKuVdd8zlctjZ2cHCwgIUCgUt92SzWYRCoZZkPYVgf38fg4OD1MXa29ujD/PY2JgoMiKpVArn5+eYmZkRzU0k5SsSZiiVSkpKr85SEgGxcrkMv9+PcDhM1R8ymQz0ev297y8SiSAej9Oxi2LiV37lV/D27dufh8Ro3Qu0zc7Zasyp0+kwOjoKj8fDm1VE5kPyiTnJ8COVSkUNs1Ao4N27d1hcXHw0lwoAdRFnZ2drfk9iXdLjubGxQWe09PT0CHqYWJaFx+PB4OBg3Y6ch4A0kyuVSiwuLlLFApVKhZ2dHVgsFlgsFuh0Orqrdnd300b1crmMy8tLpFIpdHV1UY5y9Xf8GILS+/v7vNQyHhttY5zATfZLSKq9XC5je3sbJpNJkJpeKpXiZZyEczoyMgK/3091jlwuF1UyeCx4vV6EQqE7vaUEEomEuv/r6+vIZDJ0ZMLOzg70ej3sdjt1yeuhUqlgb28PBoPhUd9LJpOh6g/kM19eXkY0GkUmk4HBYMDx8TGsViuMRiMNOeqpP+zs7EAqlaKrqwu9vb0Nhwy1ApJHAJoLSj822iohROhjfJDL5bCxsQGbzSZITQ/gx8e9vr6mzBubzQaTyYRyuQyPx4Pu7u5H7fNjWRbFYhHLy8u8qIlkWjPZcaampqBUKmmizOv11k26xONxWK1WQdzhVhAMBjE7O1uzGBJ5lf7+fkilUjidTqTTaZydnUGpVCIYDFK2EtkdifrD9PQ01Go1XVx8Ph9isZgo90quVd2T+1RoG+NkGAbJZBI+n6/psSRmJNlCoa6wyWRqOJPk6uoKPp8PKysr1N1iGAaHh4cYGxsDwzDY39+n0pat6BjVA8uyNJ4cGxtreTdTKBTo7e2liTGtVotAIIBXr16BZVmEQiHs7e3RcRKPhUQiAY/Hg+Hh4YZJPsLnHRsbw8zMDI2hT09PcXBwAPXot3kAACAASURBVIVCQemThKhPsqjd3d2QSCS0idzlcj1IUlSpVGJ2dhalUulRpr8JQVu5tXxizsvLS1xdXeHZs2c0BhQa98lksrpkB5ZlkUgkYLfba0awn52dIZ/PY3Z2FhKJBFqtll6T6AIVCgWsr68jn8+3NEKeuOh2u13Uqc8AqEQHx3HI5XI4ODiATCbD69evaaeM2KT8WCyGw8PDloY8SSQSKo9Jdi+/30+1kaamplAqlSgjS6vVoru7GyzLYmBgAJFIBEdHR5ibm8PV1RUMBgNvSdFcLkd1hJ9SBQFoM+OsbnOqB+Kera2t1ajSCd1hkskkEokExsfH6e+KxSLevn2Lrq4uSsUjs0ucTieUSmXdL5d0exB36OrqCuFwGAaDAZOTk1SMqxHIgjQ4OAibzSbovQhFOp3G0tISDAYDCoUCZDIZHYdgsVhEE2hOJBJYXl5+sAo7+eymp6fp4iKTybC7u4tyuQyZTEYXSpZlabafyNfIZDJ4vV7kcjmsrq42lVMtFos00dYxzipoNJq6X2a5XEYsFqM6rrcf9qOjI0EPdb0659HREYaGhmC1WmlGdnd3F4ODg7wabckCMT4+jrGxMSSTSUilUrhcLpRKJdjtdnR3d995KNLpNPb29rC6uvqohhmLxeDz+Woyv+SzJsOGotEoJBIJjo6OaBeKUO5uIBBAPp9/lCw2wzDUYIgnoFQqoVarsbu7C5lMVpNUIk0I5HMl/N/Dw0PI5XLMz8/fGVFI3Oan1g8C2sg4GYZBPp/HwcFBDfMim83i7du3GBkZEc31UqlUdIcm081mZmZqhLIODg4wMTHBW6v29nshyY/Z2VlkMhk6tfv09JQybCqVClwuF+bm5kRVt7uNUCgEt9vd0MUk9wTcLDDJZBKhUAgGg4Gq6tvt9oY1ZJ/PB5/P96B5pXxA+mYXFhboojg7O4t4PI5wOAybzQaPxwO5XA6z2UwbHYj6w8DAAEqlEqRSKQ4PD+msmqGhIej1eqhUKpycnDypfhDQRiQEQlTf3t6uaSfa399HX19fw7LH8fFxjYvaDMQIg8EgvF4vFhcX6RcYiUSgUql4xyhCQeRD/H4/JiYm6AgKo9H4KNdjWRbRaBRGo7HlBSCfzyMUCiEUCmF0dJTGgdX3zHEcLi8v0dvb+2i9q9XXcTqdDcOZTCZDNZLm5+fpCEKdTlejjEgy94lEAhaLBQcHB8hms/iP//gPhMNh/PM///OjvZcqtDcJAUANa+fi4gIKheJOAb4ehHJN0+k0PB4PlEolpVSRL53UFh+rtKBQKOisEJ1Oh1KphKurKxwcHGB+fp6KJYvxgHs8HhSLxQd3r6hUqpr5JfF4nN7z1NQUnWQtliL8fSBlFj4Dr4jiHnk2JBIJHeFB3NnqkhrJM5Aum+vra/z3f/83vvzlL+Of/umffq6KgQRttXMSpe13796hVCphbm6O10P6k5/85M506/tQqVTw9u1bcByH5eVlSsUrlUpwu901g3IfA5eXl4jH43docuR7uL6+htfrhUqlwvj4eMsc2pOTE+TzeVHpeLdBOnESiQTK5TIUCgVGR0epIr+YOD09RTabpRnzVlH9ORP5lYmJCVo2IV0/P/nJT7C3t4ePP/4YOp2urjsvRD+oybF1P6y2MU7g8+Zlh8MBp9PJ+wvmOw+TZGQtFgvy+TwmJydRqVTw7t07TE5OPmpdi7iXZrP5TgNyPZAmAJ/Ph2AwSKl5zbKfJMucy+UaclIfCvLcEGU6wjWWSqWIxWLwer2wWq3o6el5UGKFXMfn86G3t1f095PP5yGTyWi8XCwWEY/H8bd/+7f4yle+gr/7u7+797Wrq6tUPygej9+rH8Tj2PZ2a1mWxf/+7/+iUqkILoxrNJqmHEsiwTE6OkoFrHZ3d5FIJNDf3/+oCRlCO3M4HLwzsiQZQQY2kd7JSCSCdDoNh8NxJ04ljJlGItxigBT7SQaagBhhd3c3rFYrIpEISqUSwuEwXWDq9ag2ug7ZyRpNEn8ISFscx3Gw2+2Ix+P4kz/5E+Tzeco8qkdYEaIfJOTYarQNQ0gikWB9fb0ll6XZSINIJIKNjQ3ahU/S7ENDQ+jt7aXskqurK8TjcdFpW263uyZmEwqFQgGn0wmdToeuri6YzWZaT00kEggGg5TEYLVaH9UwyaKmVqsbSqyQHlWj0QiLxUJ7VAOBADKZDC4vL1EoFO59PSF3lMvlR9eM9Xq9iMfjcDqdePHiBf70T/8UFxcX+Pjjj+9NRDbSD3rIsdVom50TuPlC+WrJVmN/fx+jo6N1m4PJ4NWVlRXaNE00XFdXV2s+fKVSSRMdz549o2yfVmOcdDqNVCqFiYkJ0dwxUvIghpHJZHB9fY3T01P09fXB4XDw6kNsBUR+ROiovds9qoVCgcqU2mw2OJ1O5PN5miEni6Ner6eSIY8FMoh3amoK3/jGN/DlL38Z3/zmN8EwzL1NB4Aw/SChWkMEbWWcwE1sMTQ0JMggyuXyHToex3GIxWIwmUxYW1ujROZAIICrq6u6GVmbzUbFvBiGQTAYhN/vp2MhZDIZ7wclGo3i6OgI8/Pzj/pwkTGI09PTMBgMSCaTOD4+BnAjhi1EA7gRyM48MjLyYIFppVJJVdYJ6+fs7AyZTAZ9fX1Udf6xhyeR3XtwcBDf+MY3sLq6ir/4i7/g9XkJ0Q9qVWuorYyTzO8YGBgQZJy3u0zIqqxWq2takHK5HEwmExwOR8O4h3w5w8PDGBoaQjqdpgODcrkcHA4Huru773WlOY5DOBwWhb7WCOR9zszMUC1Yk8mE9fV1mvmORqM4Pj6mD3srvN1KpYI3b96gv7+/YcNAKyCsHyIM5nK5IJVKcXZ2hkAgQCd+iY2rqyvaLfNHf/RHmJ2dxbe//W3eC5kQ/aBWtIaANjNO4MbQCGeSLwiHleD09JTGXoRwsL+/D4PBIFhCkWEY+kBPT09T9QOO43B2dkaV4LVaLVV912g0j6qMB9xk/DQaDdbW1uouNMStValUMJlMlAmVSCRq5ss0S84QNXmiMvhY4DgOEokEXV1dNS57sVhEPp+nLrDdbm+o68P3WjKZDHNzc/jzP/9zDA0N4a/+6q8EeRh89IMsFgtMJlPTY+9DW5VSSGbPYDAIMs7r62v6hSWTSfT09FC6FmnvslgsoidKCNsnGAxiYmICx8fH4DgOCwsLj8qSqabjCZ2yxnEcEokEHXLU1dWFSCRC1dyrkc/n8fbtW8zNzT3YIJrd0/7+fo1h3kaxWEQ4HEapVEJfXx9OTk6o3IkQL4u0+I2MjOCb3/wmzGYzvvvd77aUV9ja2sLLly/r1i6//vWvU/2gZsfiF6HOWSqVEI1GodFoBCU0yBjzYDCI+fl5qNVqcByHTCaDcrkMo9H4qMQC4jafnJygUCggnU5jfn6e6tSKee1CoYC9vT0sLCyIUv4pFovw+XxUt2d4eBj5fB5yuRxbW1uYmppqiV/MF8Sr0Wq1DYf63H5NNBqlOrqzs7MIBoOwWCwNPxOikL+wsIC//Mu/hFKpxN///d8/6rPBE+1vnEQIqqenR1Bsc3JygmQyidnZWUilUnAch2g0ipOTE8zPzz8qgblUKmF7extTU1N0dyHG6vP5cHFxAb1ej5H/v71zD4rqPP/4d3chLBfdlcuCynW5ekEQMImNzs+JaJpExVaiNrWTtrbYptO0dTLR6cROrpNiO2basVPx0s6YZhLDorFqM1Vsx8bGGyyCXATZZZHLwsKyCwt7Y3ff3x943u7iAnvgHLJcPjOMLhw4Zy/PeZ/3uXwfuXzKXQ4dHR2IiopipRjBBkIInVXJdNIwBsNX3S8TZ5jKvpLRQert7YVYLMaqVatgtVo9vAomGBgeHo5Dhw7B6XTi6NGj/mCYgL8XITCwGenHSFOGhIQgIyODyksAI825OTk5vFb92Gw2VFVVITk52cPtY97w2NhYLF26FAMDAwgICKDdFEwzt68rH5OMt9lsrAW82CAQCGC322G322kawW63Q6lUYuHChbS4gAuYfGlcXNyUAz5CoZBKczocDvq5sFgsiIiIwMKFC9Hf34+UlBT85je/gdVqRUlJib8Y5pj43cqp1WoRHBw8YfDB6XSipqaGdsHrdDpERESgu7sbMTExvAYvANBGZbPZzOpcTPsYU0LIRH/Ha8UymUzQarVITU3lNS1jsVhQXV2NrKwsj+th9qlWqxWRkZG0CsnbPtUXXC4XVd/js1je6XSira0NOp0OZWVluHXrFkJCQnD27FmPih0/wP/dWqfTSQcBjRcQcjgcVN4/IiICw8PDdGDRE088gTVr1tA5KXx8mPV6PR48eIC8vLwp9ZjabDYaUIqLi6NliEwy3ul0orm5GSkpKbwGmICRG0BYWBgcDseEKzpzg2Hqk5kRC75sHxhX1mAwTIuBaLVaRERE4He/+x3dq3d2duLUqVO8n5sFM8M41Wo1hELhmPW1JpMJNTU1yM3NpdOugP8Fk2JiYjA8PIympiaYTCbExMQgKSmJE+FlYCR5zEwt49pl7u/vR2trK4aGhpCUlIT29naqUMcnOp0OGo3GZ7U/d1wuF7RaLW0mz87OpoY++vV2Op20kIHPIBMwcgPV6/VITU3F73//ezQ0NOBvf/sbbwOMp8jMME6NRgOHw+FV5qK/vx91dXVYtWoVxGIxlf6vr69HTk7OYy88IwAVGBiIioqKKe2bCCEwGo002c/nSuZyuaiSOZNaiouL43ToLMPAwAAaGxuRnZ095egv81lqaGhAf38/1RNmapnv3r3rIRrNF0zhRXZ2No4dO4aKigp8+umnvDY3TJGZYZxarRY2m+2xldNgMEAikVB5CWYf1NjY6FMejjnebDYjOjoaDQ0NNBE/0d2UEIL79++DEIJly5bxuu8zm81QqVRU1pK5bpFIBJvNho6ODtrdMtUPm16vR3h4OFwuF+c3G+bGGRISgsrKSqpCwEiL8olWq8WiRYvwl7/8BV9++SVKS0u/dpnLCfB/42Rma3hcwCPDYJqvGR1Zq9VKo21sS+QYSZTu7m7YbDYsX74c3d3ddDDsaFpbW+FwOCCXy3n9YJlMJtTW1mLFihVeV0kmd6vT6WjQjIn+si1G0Gg0MBqNHjo8fOBwOHD37l3ExcXRAoLa2lqEhYXRyiquMBqN6O7uRlpaGk6ePInLly/j7NmznJRQbtq0yWMytTfYNF+PYmYYZ09PD3p7e+n4Pnc3l7lWtVoNi8XC2TwLl8uF1tZW9PT0IDg4GJmZmRgeHqY1snymLxiY3KxYLPb5A2u329HV1YWenh5ERERg6dKlsNls445gAEZWTEaJj0/DdDqdUCqViI+P92j7Yqp9ent7sXLlSmi1WoSGhk5JR8loNNL5NZ9++ikuXLiAzz//nPVNazTl5eVQq9XYt2/fhK2EbJqvRzEzjLO3txft7e1ITk6GTqdDbGysRymeSqWCy+XiLa3gcDggEolQUVFBA0psJ2GzpaenB93d3VMaJMS4kc3NzbBYLIiPj6ddHe4iXF1dXbRBms8bDvM6mkymCffKvb296O7uxsDAgMfQWjautlarhVQqRWlpKRQKBf7+979zKm3p3srmDaVSiQMHDnisrosWLfJ1TMTMKEIQCoWwWq2orKykQsKMuzs4OIikpCRegzHM/nPRokU0yksIQVVVFSQSCaKjoydcmdjAlZwkM3YhKysLTqeTiiPfu3cPEokEixcvhlarhUgk4r1Hcnh4GFVVVbSNbSJGt+rp9XrU1dXR0QgBAQFjrvADAwPQarVIT0/HJ598gjNnzuDixYvTrjk7XkP1ZHqUAT8zTkaRjpHcDwoKgsvlgtlsRm1t7bTk+3p7e2G1Wh+bip2ZmYne3l60tbUhIyMDDx8+hEQimVIu1eFwICQkBKtXr+b0eYlEIlpE8NRTT8FoNKKnpwdCoRA2mw3d3d2IjIzkJa3gcDhQVVU14XwUbzCvY2JiIhITEzE0NITAwEA0NTVhYGAAkZGRiI2NpXEBZuLbqlWroFAocPr0aVy8ePFr0ZudbEP1ePiVcVosFnR2dkIkEiEoKIjeSTUaDVasWMF71U9nZyc6Ojq8rmLMOHjGVRSLxWhra4PZbKYj+EJCQnzawzESHEKh8LGbANcwwS9G13dwcBA6nY56IVqtFlFRUZwETZj2smXLlnHyXjFGlp6eTvepwEhLIBMYXLZsGS5fvozjx4/j0qVLvH9GxmKyDdXj4VfGGRwcjOjoaHz11VcARnoSg4KCOB/qOhpmL2G325GTk+PTKiaTySCTyegNpLOzE3q9HgsXLqTplrFW1MbGRggEAl6mMbvDqBcw6njAiPQH8wF2OBy0YTs0NBQZGRmTHkNgs9lw9+5d3m6izOQ0YMQNbmlpQV9fH7Zt2waLxYIjR458remSyTZUj4dfVf6Wl5fjnXfeQXBwMP7973+jo6MDWq0WLS0tMJvNvJzT5XKhoaEBer0eiYmJrN1L5kOflpaGp59+GnFxcRAKhaipqUFNTQ26u7tpFRMzGCk+Pp5TXaGx6OrqwpIlS8astgoICEBCQgLy8vKQnp4Oh8OBxsZG3Lx5k8729EXsjNljpqam8tr3CYyUDtbX1yMlJQXDw8OIiYnB6dOnUVtbi/b2dl7PPRq1Wk1Xy8k2VI+H6K233hrv5+P+kGuio6NhMBjw2muvoa6uDlarlUpZtrW1ISoqipaJcXGXZFaNsLAwTjRRBQIBdQ+ZwBEjrdjc3IympiYEBgZypuszFna7Ha2trUhISPB5FRMIBBCJRFi8eDGVYAkICMDt27cxODgIoVDotTjfbrcjMDAQ4eHhE04K5wK9Xo+4uDjcunUL7733Hi5duoSsrCzk5+dz1jHjjlKpxEcffYTy8nJYLBYIBALaRldUVASr1UqnBqxZswYfffQRLBYLrl69irffftvXVM7b3r7pV6kUADh37hw6Ozvxgx/8AF988QUUCgXq6urwf//3f9i+fTsSExPpJKu8vDy4XK5J9TcyxfMGg4GXN9UdQgju3LmD4OBgWK1WZGdnY3BwcMypalPBarWiuroaKSkpnDwvl8sFg8GAoaEhxMXFoaGhAREREYiMjITdbkd1dTXvPbPASOXUw4cPkZGRgWvXruHNN9/EpUuXPHRzZzD+n+ccC6vVisuXL0OhUKCqqgrr1q1DQUEBvvGNb0CtVqOvr4+OT/clAmk2m1FTU4NVq1bxHnI3m80QCoUQiUQe1UdMCoUQgqysLIhEIk4itkxzNx+F5YQQOn3MZDLBZrMhOjoaS5cu5XW/x7SyLV++HDU1NTh48CAuXrzIqz7vNDNzjdMdu92Oq1evorS0FLdv38batWtRUFCAjIwMxMTEoKmpCUKh0KsiOjBiLNXV1Vi5ciXvkT2mHG+8czFu4YMHD2AwGBAZGYn4+HjWdbNMOSLfQSZgJPjDjILX6XTo6elBQEAAsrOzeRnXrtVqERISgvr6erz++uu4cOEC66kAfs7sME53hoeHce3aNSgUCly/fh1r1qxBQUEBMjMzodfrkZaWBoPBgCeeeAJSqRSDg4MIDQ2F3W6fclnXRDidTlRWVmLFihU+u3yMwFlUVBQ0Gg0IIZDJZBPmUgcGBlBXV4fMzMxpCcjU1NQgOzvbYw/KKCbW1NTAYrEgMjJyyrXIVqsVGo0GGRkZuHPnDn7xi1/g/PnzrKfKeYNNHaxSqURFRQWAkbK8wsJCn/WOfGT2Gac7DocD169fh0KhwLVr15CVlYWCggLk5eVBr9dDp9MhJCQEK1asmJKKuy/09vYiLCxszFH1vuB0Oul1JyUl0SFB3tTm2traEBERwbuLbrPZoFQqJ/Q6mEnkUVFRqK6uhlgshkwmg1Qq9fn1YCRgMjIy0NLSgldffRXnzp3jzCjY1MEePnwYb7zxBn28b98+lJSUcHIdj5jdxumO0+nEjRs3oFAoUF5eTmX9//znP9Px6xKJBMnJyZwHZNzL8bjsHzQajdBqtTAajcjMzERAQABMJhMGBgamxZW1Wq0ICgqCzWZj5XW4K+WlpaWhra0NISEhEw406urqQlBQENra2rBv3z4oFApWA5LHg20dbHJyMiorK+nKOl3G6VdFCFwhEomwbt06rFu3Dq+99hq6uroQFxeH5557DikpKdi+fTsSEhIQEBCABw8e0MBGRETElFZUm81Gld65LjOUSqWQSqU079jY2IiOjg5ERETQvkm+YPbOq1evZr0dEAqFtHYWGHke3d3dUKlUVE4mJCSE3sjsdjtaWlqQlpaG+vp6FBUV4bPPPuPMMAH2dbAHDhxAUlISiouLAYD+yzez0jjd2b9/PxISEiAQCOByuXD37l0oFAp8+OGHiIuLQ0FBAdavXw+DwYDQ0FCYTCYAI1UobEbVaTQaLF26dNzhN1zAdEeEhYVh/fr1VGiM6c+UyWRUPpMLLBYL1d7hYp8ukUg88qGMEoNIJEJGRgbu3buH1NRU3L9/H3v37sUnn3yCjIyMKZ/XHbZldUVFRTAajXS1zM/P511mBZilbq0vEEJQW1uL0tJS/OMf/4BMJkNBQQGeffZZWCwW9Pb20o6IwMDAMVM0TIURMyGZ76qf9vb2MWVcmEbsyMhIWCwW2Gw2yGSySbvuZrMZwcHBsNvtvM58AUbcZqPRiAcPHmD//v0wm804cuQItm/fzvlrqlAoUFJS8phbe/XqVa8rp/ue8/jx4yguLoZKpeLykuaOW+sLAoEAmZmZyMzMxNtvv4379+9DoVBgz549kEgkKCgogFwuh9PpxMOHDxEUFOR1DPzQ0BBCQ0Pp6swnra2t6OvrG3N1Dg0NpbNgAgMDodPpcO/ePSQmJlK311f3t7+/Hw0NDbwPYwJGotStra1IS0uD2WyGWCzG7t278fHHH0MikeDZZ5/l9Hxs6mCZEQoMRUVFUKlUU2oF85U5a5zuMN0Nhw4dwptvvonm5maUlZXh5ZdfhlgsRkFBATZt2oTAwEBqIJGRkTAajdPSLcN4N0FBQcjKyvJpXywWixEfH091Yfv6+qDRaDA8PIy0tDSaA/Z2QzGbzWhoaEBWVta0GObdu3eRkJCAhw8f4rvf/S5OnjyJNWvW8HZONkOIwsPDUV5e/lidLN+GCcxht9YXmL1kWVkZzp8/D6FQiK1btyI3Nxf379/Hhg0bqNvHVduVt2t48OABQkJCOEm8MxpNer0eGo2GTmNjcrFM6xsfxQTe6O7upkrzO3fuxLFjx/D000/zfl42Q4gUCgVdaY1GI/Lz87k2zrmTSuEDZix9SUkJjh49itWrV2Pz5s148cUXERAQAJ1OB7lcTnOoXARPGHEzoVDIy36WyaWKxWKq7Dc4OIg1a9bwvmI6HA6o1WqkpqZCq9XipZdewh//+EesX7+e1/P6KfN7zqkgEAgQGxtLOw6WLFmCs2fP4le/+hUGBwexZcsWhIeH09WWEILU1NRJi1YxbWYymQzh4eG87GdFIhFkMhmAkXTJ0NAQneW5cOFCmM3mKaeXvMGISy9ZsgTd3d3YuXMnjhw5MlcNc0zmV04O6Onpwblz51BWVoa+vj688MIL2Lp1K5KTk6HX6+mw2sWLF/u0ojJDfhYvXkyNh0+8jWIwm81ob2+HXq9HQkICvQ4upE16enqoasKOHTtQXFw85d5HBrbylAqFwuNxYWEhJ9fBknm3djro6+vD+fPnUVZWhq6uLmzevBkFBQVYtGgRFixYALvdDpPJBJlM5rUOlsnFRkVF8T6GARgxlJaWljFHMTACa4ODg2hoaKCBpsnMOXE6nWhpaaE3rW9/+9t477338M1vfpOLpwKAfVmeXC5HYWEhjEYjNm7cSH93mpk3zummv78fFy5cQFlZGTQaDTZt2oRt27YhOjoaPT09iImJgUQigcPhQFhYGJxOJ4RCIYaGhqZFC8dkMtF0ia9FC0NDQ3C5XHA6nVCpVFSuZaI9qsvlQnV1NaKiohASEoIdO3bg0KFD2LJlCxdPBQC7sjyj0YikpCRfpSv5Zn7POd1IJBLs2bMHe/bsgclkwqVLl/CHP/wBTU1N2LhxI7Zv347g4GA6vMjlcmHZsmWsBgdPlv7+fixcuBC5ubmsSg3dO2yWL1+Onp4e6PV6SKVS6HQ6yGQyr7lUJv0UFhaGHTt24ODBg5waJsCuLK+iogJyuRwKhQJSqRRKpZKPbpMpwbtxstkDTEHO3u9ZsGABdu/ejd27d8NsNuOLL77AsWPHcO/ePTz11FOor6/HyZMn4XK5cPPmTZri4KMFrKurC21tbVi9evWU9pDBwcE0j8rsIe/fvw+JRILExERYLBYEBwdDo9EgKSkJg4ODKCwsxP79+/Gtb32Lq6dDYVOWp1aroVQqaSleXl4ecnNzua78mRK8u7Vs9gBTkLOfsTQ2NuKFF17AypUroVarqcpDWloalQfV6XSIjo5m1XI1FgMDA2hqakJ2djZv4/AIIbBarWhsbITBYIDdbkdwcDDef/99vPrqq/jOd77Dy3nZlOWVl5dj3759HsbITBT4GlbP6XdrlUqlR+BAKpWivLx8ysfONs6cOYO8vDwPlYdbt25RlYcVK1ZAq9VicHAQUVFRMJvNXvs6J8JgMEAqlfLSNeOOQCBAcHAwHVvY1dWFV155BSKRCDdu3MCLL77IyzhDNmV53gzQ37w0XqUxx9sDTOXY2UR6ejry8vIAjCgKPv/88zh16hSqqqqwa9cuXLx4EVu3bsXRo0fR0NBAFdtv3bqFjo4OuFwumhMdj46ODrS0tPAy7m80hBC0tLRg0aJFWLx4MT744AO88847UKlU2LZtG29iYL6U5TFSlswAX+ax0WiEXC6fO3tONnsAPuTsZzKBgYHIz89Hfn6+h8rDoUOHqMrDypUrMTg4iLq6OixYsABLly71GkxiRuMxQmJ8QghBfX09xGIx7HY7Xn75ZezcuRM//OEPIRAIOMtnjsWJEydw+PBhi+1jMAAAB5BJREFUGrc4ceIE/dmBAwc8yvJKS0vxwQcfIDk5GSqVyu+2ULwaJxuJej7k7GcLAQEB2LBhAzZs2OCh8vDuu+9i2bJl2L59O73jm0wmaDQaOmDXaDQiPDwc2dnZvEqzMBiNRgQFBSE2NhZ79uzB1q1bUVRUxHvHDkNOTg5dQUffCEYbn1wun7bG6cnAq3Gy3QNwLWc/G3FXeXC5XKioqEBpaSl++9vfUpWHZ555BgMDA2htbUVAQAAtzOfTOAkhePjwIWJjYxEaGopXXnkF+fn5+NnPfjZthjnb4NU42bTm+CpnPxnVNKPRiDt37qC4uNiv9hRTRSgU4sknn8STTz6J4uJiVFdXo7S0FEeOHIFYLIZEIsGxY8dgMplQUVEBmUyGxMREOJ1OTvWNCCFobGyEUCiE0+nE3r178cwzz+CXv/wlJ4Y52RQbD1o/0wszf3KMrylTWVlJiouLyZUrV0hxcTExGAz0Z4WFhaSkpMSnYxlycnLo/w0GAyksLPR6XoPB4PG3r1y5QuRyORdPye/58MMPyebNm8mvf/1rkpubS55//nly7NgxotFoSG9vL/nXv/5Frl+/TlpaWsjQ0NCUvzo6OohSqSRGo5Hs3LmTvPvuu8TlcnH2fHx9z925cuUKGfl4zwi82t+MKt9jU56lVCrx0ksv0TyW0Wikx/pbyJxrmpubkZiYiICAANp2plAocPHiRarysGnTJoSGhiI4OBh1dXWIjo6GTCZj1cNJCEF7eztVXv/JT36ClJQUvPXWW5y5spOZGG00GqFWq7Fx40Z/Kc+bCK8vll9NGZsINumWnJwcjze0oqKCKtjNdlJSUmiBgbvKw82bN/GnP/0JJpMJ3//+9/G9730Pp0+fRlRUFBwOB9ra2qgsi9VqHfcchBA0NzfDZDKBEIKf//zniI+P59Qwgcml2EZLi8xUZpRxso3euu8vS0pKPMLqcxGBQIDU1FQcPHgQ//3vf3Hq1Ck4HA786Ec/wk9/+lP885//RGdnJwCgtrYWDQ0NtNpnNCaTCcPDw0hPT8f+/fsRGRmJ999/n/PgD9v33JukyExlRhnnZNMtx48fx65du76uXj2/RCAQICkpCa+//jq+/PJLfPzxxwgMDERRURH27t2Lr776isqV1NfX4/bt23T+ZWdnJ0JCQpCRkYE33ngDoaGhOHz4MC/RYDbvuXuAcTYw4/acP/7xjz167ibafzAlgGPdTedsJHAMCCHo7u7G2bNncfbsWarysGXLFkRFRdEqpf/85z/o6emBWCzG0aNHeUvTsHnP3bV+gP+9R/n5+f4epffubowVKSIcRWu5xj1yp1KpPCJ3KpXKI8JbWVlJKisr6ePS0tJx/94sjQROCZ1OR0pKSshzzz1H0tPTybp168iNGzfIrl27iFwuJ7m5ueTzzz/n9RrYvOfuzKD3yKv9zTjj9DU1o1KpCEZWfvo1OpVSWVlJ8vPzPb4nlUrHPb/BYCCVlZUTHjfb+Otf/0q2bdtGTpw4QdauXUvWrl1LHA4HGRgYIJ2dnbyem006jpCR96i4uJgAIMXFxUSlUvF6fRww81MpXKNQKHDmzBmPsq7k5GSUlpaOGe1TKBQoLCyc0J2ebfT29kIikdDiBULIfOUPd8wrIYxmLkcC2cIMImKYN0z+mVHRWq6Zy5HAefyfOb1ysim2VyqV6Ovr85hwfPz48ZkQCfRL5mukfWCszSjx04AQ18yBSKBfMl8j7YFX+5vTbi3wv+bc8vJyKBSKx5pzP/vsM4/jjUYjDh8+DGBE91StVnv8XK1W0793+PDhx9zm0SgUCo+vuQAbSRq1Wu3Rc5mXl+ehaDCrGctqyRxZObmGTd60uLiY5l4NBoPH785mSktLH3td5HK5R07aHfdUyJUrV2ZjGsur/c3pPSfXsFkRjEYjPvjgA5qOkUqlX5fa+LQzXyPtG3PereUSNh0U7qLGjAs82kWerczXSPvG/MrJIbNN1JgvJiNJU15eDrlcPqfyzPMrJ4ewWREYGUZmpZVKpVCr1XNi9WQjYQn8b7vAHDNXAmfzKyeHzDZRYz7xVcJSrVYjNzfX43eZyWCznXnj5BA2gmbuosbMv2OJGk8mYQ+MBJ38bTgPg68SlnK5fKRDYw4yUeH7PCwRCAQ5APIBKAHkADhOCDE++lkpgCuEkOOPHssB7AOgApAMoIQQ8phfKxAIKgkhuY/+LwVwghDy0hjnf4MQctjtcQkhZB+Xz3Ge6WHeOP2cR8ZeTAjZ5PY9AyHE65xAgUCgApDrdkOYN84Zyrxb6//IAYwuh+kTCAQ5hBBvKlfFAFoEAsGBR48PeDmGFx55AoXw4jVM5di5yrxx+j+s5rsTQo4/cn2Z1bIcjxs3X5S6ud8VAE4A8Op+szx2TjKfSvF/+gCMjv6MabDMnvPRB78EwJWxjuWSR+43DVU/WgW9JiXZHDuXmTdO/0cNL8bozaUVCARMIIo55jgAxSNjeAyBQDCh4QoEArlAIHhDIBDkP/p3rHzPmO73FI+ds8y7tX4OIUTprjrwaK9WPupx36PVpw8jK1D56L/h/viREcvh22rlq/vJxv1m5arPVf4fbHJP5PKtXUoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def bottom(x, on_boundary):\n",
    "    return near(x[2], 0.)\n",
    "bc = DirichletBC(W, Constant((0, 0, 0, 0, 0, 0)), bottom)\n",
    "\n",
    "u = Function(W)\n",
    "solve(k_form == l_form, u, bc)\n",
    "plot(u[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matplotlib functions cannot plot functions defined on a 1D mesh embedded in a 3D space so that we output the solution fields to XDMF format for vizualisation with Paraview for instance. We also export the bending moments by projecting them on a suitable function space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ffile = XDMFFile(\"shell-results.xdmf\")\n",
    "ffile.parameters[\"functions_share_mesh\"] = True\n",
    "v = u.sub(0, True)\n",
    "v.rename(\"Displacement\", v.name())\n",
    "ffile.write(v, 0.)\n",
    "theta = u.sub(1, True)\n",
    "theta.rename(\"Rotation\", theta.name())\n",
    "ffile.write(theta, 0.)\n",
    "\n",
    "V1 = VectorFunctionSpace(mesh, \"CG\", 1, dim=2)\n",
    "M = Function(V1, name=\"Bending moments (M1,M2)\")\n",
    "Sig = generalized_stresses(u)\n",
    "M.assign(project(as_vector([Sig[4], Sig[5]]), V1))\n",
    "ffile.write(M, 0)\n",
    "\n",
    "ffile.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bending moment amplitude over the structure would look like this:"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. figure:: shell_moments.png\n",
    "   :scale: 50%\n",
    "   :align: center"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Format de la Cellule Texte Brut",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
